Analiza zbioru danych Bodyfat

Słowem wstępu

Dane pochodzą ze strony internatowej Kaggle.com ze zbioru danych Bodyfat.

Dane wymieniają szacunkowe wartości procentowe tkanki tłuszczowej mierzone pod wodą, wagę i pomiary różnych obwodów ciała 252 mężczyzn.

Rows: 252 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (15): Density, BodyFat, Age, Weight, Height, Neck, Chest, Abdomen, Hip, ...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Density	BodyFat	Age	Weight	Height	Neck	Chest	Abdomen	Hip	Thigh	Knee	Ankle	Biceps	Forearm	Wrist
1.07	12.3	23	154.25	67.75	36.2	93.1	85.2	94.5	59.0	37.3	21.9	32.0	27.4	17.1
1.09	6.1	22	173.25	72.25	38.5	93.6	83.0	98.7	58.7	37.3	23.4	30.5	28.9	18.2
1.04	25.3	22	154.00	66.25	34.0	95.8	87.9	99.2	59.6	38.9	24.0	28.8	25.2	16.6
1.08	10.4	26	184.75	72.25	37.4	101.8	86.4	101.2	60.1	37.3	22.8	32.4	29.4	18.2
1.03	28.7	24	184.25	71.25	34.4	97.3	100.0	101.9	63.2	42.2	24.0	32.2	27.7	17.7
1.05	20.9	24	210.25	74.75	39.0	104.5	94.4	107.8	66.0	42.0	25.6	35.7	30.6	18.8

Objaśnienie zmiennych zbioru danych:

• Density - Gęstość wyznaczona na podstawie ważenia pod wodą

• BodyFat - Procent tkanki tłuszczowej z równania Siri (1956)

• Age - Wiek (lata)

• Weight - Waga (w funtach)

• Height - Wysokość (cale)

• Neck - Obwód szyi (cm)

• Chest - Obwód klatki piersiowej (cm)

• Abdomen - Obwód brzucha (cm)

• Hip - Obwód bioder (cm)

• Thigh - Obwód uda (cm)

• Knee - Obwód kolana (cm)

• Ankle - Obwód kostki (cm)

• Biceps - Obwód bicepsa (cm)

• Forearm - Obwód przedramienia (cm)

• Wrist - Obwód nadgarstka (cm)

Równanie Siri:

\hspace{0.8cm} Fat(\%) = \frac{495}{density}-450

Przekształcimy wagę z futów na kilogramy i wysokość w calach na cm według standardu środkowo-europejskigo.

Code

dane["Weight2"] <- round(dane$Weight*0.435,2)

dane["Height2"] <- round(dane$Height*2.54,2)

Dodamy kolumnę ID dla łatwiejszej obróbki danych.

Code

ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)

Również dodamy nową kolumnę BMI - to wskaźnik masy ciała, który każdy jest w stanie wyznaczyć i zinterpretować w której kategorii znajduje się jego wynik: niedowaga, waga prawidłowa, nadwaga i otyłość. Uważamy, że nowa zmienna BMI będzie dobrym dopełnieniem naszego zbioru danych.

Code

dane["BMI"] <- round(dane$Weight2/(dane$Height2/100)^2)

Dodamy kolumnę category dla zmiennej BMI, która ma 3 stany:

• za niska
• w normie
• za wysoka

Code

dane["Category"] <- as.factor(ifelse(dane$BMI< 20.1, '1',
                              ifelse(dane$BMI< 25.9, '2','3')))

Dodamy też kategorie dla zmiennej BodyFat.

Code

dane["Category2"] <- as.factor(ifelse(dane$BodyFat< 11, '1',
                              ifelse(dane$BodyFat< 24, '2','3')))

Sprawdzimy czy mamy braki danych w poszczególnych kolumnach:

Code

apply(dane, 2, function(x) sum(is.na(x)))

       ID   Density   BodyFat       Age    Weight    Height      Neck     Chest 
        0         0         0         0         0         0         0         0 
  Abdomen       Hip     Thigh      Knee     Ankle    Biceps   Forearm     Wrist 
        0         0         0         0         0         0         0         0 
  Weight2   Height2       BMI  Category Category2 
        0         0         0         0         0 

Nie mamy braków danych w zbiorze. Kolejnym krokiem będzie wyznaczenie podstawowych statystyk (bez kolumny ID):

Code

summary(dane[,-1])

    Density         BodyFat           Age            Weight     
 Min.   :0.995   Min.   : 0.00   Min.   :22.00   Min.   :118.5  
 1st Qu.:1.041   1st Qu.:12.47   1st Qu.:35.75   1st Qu.:159.0  
 Median :1.055   Median :19.20   Median :43.00   Median :176.5  
 Mean   :1.056   Mean   :19.15   Mean   :44.88   Mean   :178.9  
 3rd Qu.:1.070   3rd Qu.:25.30   3rd Qu.:54.00   3rd Qu.:197.0  
 Max.   :1.109   Max.   :47.50   Max.   :81.00   Max.   :363.1  
     Height           Neck           Chest           Abdomen      
 Min.   :29.50   Min.   :31.10   Min.   : 79.30   Min.   : 69.40  
 1st Qu.:68.25   1st Qu.:36.40   1st Qu.: 94.35   1st Qu.: 84.58  
 Median :70.00   Median :38.00   Median : 99.65   Median : 90.95  
 Mean   :70.15   Mean   :37.99   Mean   :100.82   Mean   : 92.56  
 3rd Qu.:72.25   3rd Qu.:39.42   3rd Qu.:105.38   3rd Qu.: 99.33  
 Max.   :77.75   Max.   :51.20   Max.   :136.20   Max.   :148.10  
      Hip            Thigh            Knee           Ankle          Biceps     
 Min.   : 85.0   Min.   :47.20   Min.   :33.00   Min.   :19.1   Min.   :24.80  
 1st Qu.: 95.5   1st Qu.:56.00   1st Qu.:36.98   1st Qu.:22.0   1st Qu.:30.20  
 Median : 99.3   Median :59.00   Median :38.50   Median :22.8   Median :32.05  
 Mean   : 99.9   Mean   :59.41   Mean   :38.59   Mean   :23.1   Mean   :32.27  
 3rd Qu.:103.5   3rd Qu.:62.35   3rd Qu.:39.92   3rd Qu.:24.0   3rd Qu.:34.33  
 Max.   :147.7   Max.   :87.30   Max.   :49.10   Max.   :33.9   Max.   :45.00  
    Forearm          Wrist          Weight2          Height2      
 Min.   :21.00   Min.   :15.80   Min.   : 51.55   Min.   : 74.93  
 1st Qu.:27.30   1st Qu.:17.60   1st Qu.: 69.16   1st Qu.:173.35  
 Median :28.70   Median :18.30   Median : 76.78   Median :177.80  
 Mean   :28.66   Mean   :18.23   Mean   : 77.83   Mean   :178.18  
 3rd Qu.:30.00   3rd Qu.:18.80   3rd Qu.: 85.69   3rd Qu.:183.52  
 Max.   :34.90   Max.   :21.40   Max.   :157.97   Max.   :197.49  
      BMI         Category Category2
 Min.   : 17.00   1: 29    1: 48    
 1st Qu.: 22.00   2:145    2:128    
 Median : 24.00   3: 78    3: 76    
 Mean   : 24.87                     
 3rd Qu.: 26.00                     
 Max.   :159.00                     

Zauważyłyśmy, że wartość minimalna zmiennej BodyFat jest równa 0. Przyjrzymy się tej obserwacji.

Code

dane %>% 
  filter(BodyFat == 0)

   ID Density BodyFat Age Weight Height Neck Chest Abdomen Hip Thigh Knee Ankle
1 182  1.1089       0  40  118.5     68 33.8  79.3    69.4  85  47.2 33.5  20.2
  Biceps Forearm Wrist Weight2 Height2 BMI Category Category2
1   27.7    24.6  16.5   51.55  172.72  17        1         1

Przypuszczamy, że to może być błąd we wprowadzeniu danych.

Code

495/(dane[182,"Density"])-450

[1] -3.611687

Obliczając poziom tłuszczu z równania Siri uzyskałyśmy wartość ujemną, co oznacza, że wzór Siri nie jest optymalny dla każdego przypadku. Ze względu na możliwą do uzyskania wysokość i wagę zostawiamy obserwację w zbiorze danych ze zmienionym poziomem tłuszczu.

Code

dane[182,2] <- 495/(dane[182,"Density"])-450

Przedstawimy podstawowe statystyki danych w bardziej zgrabnej formie:

Code

st_op <- round(apply(dane[,c(1,2,3,8,16,17,18)],2,summary),2)
st_op <- rbind(st_op,St.dev=apply(dane[,c(1,2,3,8,16,17)],2,sd))
colnames(st_op) <- c("Density", "BodyFat","Age" ,"Abdomen", "Weight", "Height", "BMI")
rownames(st_op) <- c('minimum','kwantyl dolny','mediana','średnia','kwantyl górny','maksimum','odchylenie standardowe')
as.data.frame(round(st_op,2)) %>%
kable(caption="Statystyki zbioru")

Statystyki zbioru

	Density	BodyFat	Age	Abdomen	Weight	Height	BMI
minimum	1.00	-3.61	0.00	79.30	15.80	51.55	74.93
kwantyl dolny	63.75	1.04	12.47	94.35	17.60	69.16	173.35
mediana	126.50	1.05	19.20	99.65	18.30	76.78	177.80
średnia	126.50	1.04	19.15	100.82	18.23	77.83	178.18
kwantyl górny	189.25	1.07	25.30	105.38	18.80	85.69	183.52
maksimum	252.00	1.10	47.50	136.20	21.40	157.97	197.49
odchylenie standardowe	72.89	0.29	8.37	8.43	0.93	12.78	72.89

Przedstawimy również graficzną prezentację rozkładów :

Code

ggplot(dane, aes(x=Weight2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(50, 160, 10)) +
  labs(title = "Histogram wagi [kg]", x = "Waga [kg]", y = "Liczebność") + theme_bw()

Code

ggplot(dane, aes(x=Height2)) + geom_histogram(color="gray47", fill="navajowhite2") + scale_x_continuous(breaks=seq(74, 204, 10)) +
  labs(title = "Histogram wzrostu [cm]", x = "Wzrost [cm]", y = "Liczebność") + theme_bw()

Na wykresach widzimy, że niektóre obserwacje znacznie wychodzą poza tendencję, w dalszym ciągu zobaczymy jak to wpłynie na wyniki.

Ze statystyk opisowych odczytałyśmy wartość maksymalną BMI, która wyniosła 159. Wyświetlimy obserwację posiadającą niemożliwy poziom wskaźnika masy ciała:

Code

dane[which(dane$BMI == 159), c("ID","Weight2", "Height2", "Age")]

   ID Weight2 Height2 Age
42 42   89.17   74.93  44

Jest ona najprawdopodobniej pomyłką we wprowadzaniu danych, przejrzymy się tej obserwacji w dalszych analizach.

Na koniec przedstawimy wykresy, który ilustrują podział wskaźnika masy ciała i poziomu tłuszczu na kategorie: za niska, w normie, za wysoka.

Code

dane[,-1] %>% 
   ggplot(aes(x = Category))+
   geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
  scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
   theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
   labs(title = "Rozkład BMI", x = "kategoria wagi", y = "liczebność")+
   scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
   guides(fill=guide_legend(title="kategoria"))+
   scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
   theme_minimal()

Najbardziej liczną jest grupa z prawidłową masą ciała. Tylko 29 osób z badanych 252 mieści się w kategorii BMI poniżej normy.

Code

dane[,-1] %>% 
   ggplot(aes(x = Category2))+
   geom_bar(color = "gray70", fill = c("wheat1", "wheat2", "wheat3"))+
  scale_fill_manual("legend", values = c("za niska" = "black", "w normie" = "orange", "za wysoka" = "blue"))+
   theme(axis.text.x = c("niska", "w normie", "za wysoka"))+
   labs(title = "Rozkład BodyFat", x = "kategoria poziomu tłuszczu", y = "liczebność")+
   scale_x_discrete(labels = c("za niska","w normie","za wysoka"))+
   guides(fill=guide_legend(title="kategoria"))+
   scale_fill_discrete(breaks=c("1", "2", "3"), labels=c("za niska", "w normie", "za wysoka"))+
   theme_minimal()

Rozkład poziomu tłuszczu w grupie badanych mężczyzn zachowuje się podobnie jak rozkład BMI.

Code

dane <- dane[,-21]

Hipotezy

1. Czy poszczególne obwody: szyja, klatka piersiowa, brzuch, biodra, uda, kolana, kostka, biceps, przedramię i nadgarstek mają wpływ na BMI?

2. Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?

Funkcje

Code

gauss_markov <- function(mod) { 
  tabelka <- NULL
  wiersz <- NULL
  bp <- bptest(mod)
  gq <- gqtest(mod)
  dw <- dwtest(mod)
  bg <- bgtest(mod, order.by = ~fitted(mod), order = 3)
  sh <- shapiro.test(mod$residuals)
  ks <- ks.test(mod$residuals, "pnorm")
  lil <- nortest::lillie.test(mod$residuals)
  res <- resettest(mod, type = "regressor")
  rai <- raintest(mod)
  har <- harvtest(mod) 
  w <- c(bp$p.value, gq$p.value, dw$p.value, bg$p.value, sh$p.value, ks$p.value, lil$p.value, res$p.value, rai$p.value, har$p.value)
  h0 <- c("jednorodność wariancji", "jednorodność wariancji", "brak autokorelacji reszt","brak autokorelacji do rzędu 3","rozklad reszt jest normalny", "rozklad reszt  jest normalny", "rozklad reszt  jest normalny", "zależność jest liniowa", "zależność jest liniowa", "zależność jest liniowa")
  h1 <- c("brak jednorodności wariancji", "brak jednorodności wariancji", "występuje autokorelacja reszt", "występuje autokorelacja rzędu 3", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "rozklad reszt nie jest normalny", "brak zależności liniowej", "brak zależności liniowej", "brak zależności liniowej")
  n <- c("Test Breuscha-Pagana", "Test Goldfelda-Quandta", "Test Durbina-Watsona", "Test Breuscha-Godfreya", "Test Shapiro-Wilka", "Test Kolmogorova-Smirnova", "Test Lillieforsa", "Test RESET Ramseya", "Test Rainbow Uttsa", "Test Harveya-Colliera")
  for (i in 1:length(w)) {
    wiersz <- c(n[i], round(w[i], 3), 
                if (w[i]<0.05) {
                  h1[i]} else {h0[i]} )
    tabelka <- rbind(tabelka, wiersz)
  }

  colnames(tabelka) <- c("Nazwa testu", "P-wartość", "Wniosek")
  rownames(tabelka) <- c(1:length(n))
  tabelka %>% 
    kable(caption="Założenia Gaussa Markowa")
}

Code

odstajace <- function(mod, dane, stopien) { 
h <- hatvalues(mod)[hatvalues(mod) > 2*3/nrow(dane)] %>% as.data.frame()
h <- rownames(h)
rsr <- rstandard(mod)[abs(rstandard(mod)) > 2] %>% as.data.frame()
rsr <- rownames(rsr)
rst <- rstudent(mod)[abs(rstudent(mod)) > 3] %>% as.data.frame()
rst <- rownames(rst)
c <- cooks.distance(mod)[cooks.distance(mod) > 4/nrow(dane)] %>% as.data.frame()
c <- rownames(c)
w <- c(h, rsr, rst, c)
b <- unique(w)
counts <- NULL
for (i in 1:length(b)) {
  counts <- c(counts, sum(w==b[i]))
}
counts <- cbind(b, counts)
a <- c("obserwacja", "ile spełnia kryteriów")
counts <- counts[which(counts[,2]>=stopien),] 
colnames(counts) <- a
counts %>%  kable(caption="Obserwacje odstające")
}

Code

odstajace_col <- function(dat){
  a <- c()
for(i in c(7:16, 19)) {
  a <- c(a, dat %>% identify_outliers(colnames(dat)[i]) %>% select(ID))
}
names(a) <- colnames(dat)[c(7:16, 19)]
print(a)
a <- plyr::ldply(a, data.frame)
repetitions <- names(which(table(a[,2])>=5))
print("Obserwacje występujące co najmniej 5 razy: ")
print(repetitions)
}

Hipoteza 1

BMI to inaczej wskaźnik masy ciała - z angielskiego body mass index - jest to współczynnik powstały przez podzielenie masy ciała podanej w kilogramach przez kwadrat wysokości podanej w metrach. Klasyfikacja wskaźnika BMI została opracowana wyłącznie dla dorosłych.

BMI = \frac{masa[kg]}{wysokość^2[m]}

Code

knitr::include_graphics("bmi.jpg")

Model pierwszy

\hat{BMI} = \beta_{0} + \beta_{1} \cdot Neck + \beta_{2} \cdot Chest + \beta_{3} \cdot Abdomen + \beta_{4} \cdot Hip + \beta_{5} \cdot Thigh + \beta_{6} \cdot Knee + \beta_{7} \cdot Ankle + \beta_{8} \cdot Biceps + \beta_{9} \cdot Forearm + \beta{10} \cdot Wrist

Code

model <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane)
summary(model)

Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Wrist, data = dane)

Residuals:
    Min      1Q  Median      3Q     Max 
 -6.252  -2.111  -0.352   1.163 120.578 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -11.04020   12.45915  -0.886   0.3764  
Neck         -0.56573    0.41818  -1.353   0.1774  
Chest         0.11726    0.17128   0.685   0.4943  
Abdomen       0.06517    0.14352   0.454   0.6501  
Hip           0.48499    0.22593   2.147   0.0328 *
Thigh         0.26534    0.25046   1.059   0.2905  
Knee          0.03295    0.42583   0.077   0.9384  
Ankle        -0.13217    0.40735  -0.324   0.7459  
Biceps       -0.10045    0.31902  -0.315   0.7531  
Forearm       0.03633    0.37193   0.098   0.9223  
Wrist        -1.13466    0.94514  -1.201   0.2311  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.163 on 241 degrees of freedom
Multiple R-squared:  0.241, Adjusted R-squared:  0.2095 
F-statistic: 7.653 on 10 and 241 DF,  p-value: 1.323e-10

Code

summary(dane$BMI-model$fitted.values)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -6.252  -2.111  -0.352   0.000   1.163 120.578 

Mediana jest mniejsza od średniej to rozkład reszt jest prawostronnie asymetryczny.

Wykresy diagnostyczne:

Code

autoplot(model)

Założenia Gaussa Markowa

Code

gauss_markov(model)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.039	brak jednorodności wariancji
Test Goldfelda-Quandta	1	jednorodność wariancji
Test Durbina-Watsona	0.769	brak autokorelacji reszt
Test Breuscha-Godfreya	0.794	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0	rozklad reszt nie jest normalny
Test Kolmogorova-Smirnova	0	rozklad reszt nie jest normalny
Test Lillieforsa	0	rozklad reszt nie jest normalny
Test RESET Ramseya	0.034	brak zależności liniowej
Test Rainbow Uttsa	0	brak zależności liniowej
Test Harveya-Colliera	0.128	zależność jest liniowa

Budując model pełny otrzymałyśmy skorygowany R^2 równy 0.21, przyjrzałyśmy się wykresom i zauważyłyśmy, że obserwacja 42 jest nietypowa. Sprawdzimy, czy tak jest na prawdę.

Code

odstajace(model, dane = dane, 2)

Obserwacje odstające

obserwacja	ile spełnia kryteriów
42	4
86	2

Z powyższego zestawienia widać, że obserwacja 42 może wpływać na to, że nasz model nie spełnia założeń. Sprawdźmy czy w wyniku usuwania jej ze zbioru, nasz model się poprawi.

Model drugi

Code

dane2 <- dane[-42,]
model2 <- lm(BMI~Neck+Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Wrist, data=dane2)
summary(model2)

Call:
lm(formula = BMI ~ Neck + Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Wrist, data = dane2)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0303 -0.6120 -0.0121  0.6024  4.0308 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -12.78573    1.66370  -7.685 3.90e-13 ***
Neck          0.07917    0.05612   1.411 0.159594    
Chest         0.13840    0.02287   6.051 5.49e-09 ***
Abdomen       0.11926    0.01917   6.222 2.17e-09 ***
Hip           0.08192    0.03037   2.697 0.007485 ** 
Thigh         0.12382    0.03347   3.700 0.000268 ***
Knee         -0.24602    0.05691  -4.323 2.26e-05 ***
Ankle         0.07697    0.05442   1.414 0.158566    
Biceps        0.06284    0.04262   1.474 0.141674    
Forearm      -0.02033    0.04967  -0.409 0.682587    
Wrist        -0.00730    0.12658  -0.058 0.954061    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.09 on 240 degrees of freedom
Multiple R-squared:  0.9078,    Adjusted R-squared:  0.904 
F-statistic: 236.4 on 10 and 240 DF,  p-value: < 2.2e-16

Code

summary(dane2$BMI-model2$fitted.values)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.03026 -0.61203 -0.01207  0.00000  0.60240  4.03083 

Mediana jest bardzo bliska 0, co wskazuję, na brak asymetrii.

Sprawdzimy jak wyglądają wykresy diagnostyczne:

Code

autoplot(model2)

Z wykresów diagnostycznych raczej spodziewamy się niejednorodności wariancji i braku normalności rozkładu reszt. Sprawdzamy to za pomocą testów.

Założenia Gaussa Markowa

Code

gauss_markov(model2)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.001	brak jednorodności wariancji
Test Goldfelda-Quandta	0.753	jednorodność wariancji
Test Durbina-Watsona	0.135	brak autokorelacji reszt
Test Breuscha-Godfreya	0.575	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.024	rozklad reszt nie jest normalny
Test Kolmogorova-Smirnova	0.814	rozklad reszt jest normalny
Test Lillieforsa	0.15	rozklad reszt jest normalny
Test RESET Ramseya	0	brak zależności liniowej
Test Rainbow Uttsa	0.094	zależność jest liniowa
Test Harveya-Colliera	0.811	zależność jest liniowa

Po sprawdzeniu założeń dochodzimy do wniosku, że jednak większa część z nich na razie nie jest spełniona. Sprawdzimy obserwacje odstające dla poszczególnych zmiennych:

Code

odstajace_col(dane)

$Neck
[1]  39  45 106

$Chest
[1] 39 41

$Abdomen
[1]  39  41 216

$Hip
[1] 35 39 41

$Thigh
[1]  39  41 152 169

$Knee
[1]  39 192 244

$Ankle
[1] 31 39 86

$Biceps
[1] 39

$Forearm
[1]  45 159 175 206 226

$Wrist
[1]  39  41 226 252

$BMI
[1]  39  41  42 216

[1] "Obserwacje występujące co najmniej 5 razy: "
[1] "39" "41"

Zauważyłyśmy, że obserwacje 39, 41 są odstające prawie dla każdej zmiennej. Również obserwacja 39 jest w każdej statystyce rozpatrywanej w funkcji odstajace(). Sprawdzimy jak model będzie wygłądać bez tych obserwacji:

Code

dane3 <- dane[-c(39,41,42),]

Po wyczyszczeniu danych rozdzielamy nasz zbior na testowy i uczący.

Code

ID <- c(1:length(dane$Density))
dane <- cbind("ID" = ID, dane)
set.seed(221)
wektor <- sample(1:249, 180, replace = FALSE)
uczacy <- dane3[wektor,]
testowy <- dane3[-wektor,]

Dalsze analizy będziemy przeprowadzać na zbiorze uczącym.

Model trzeci

Code

model3 <- lm(BMI~Chest+Neck+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm, data=uczacy)
summary(model3)

Call:
lm(formula = BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6920 -0.6510  0.0483  0.6423  3.4233 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.163994   1.854714  -4.941 1.85e-06 ***
Chest        0.139161   0.025767   5.401 2.21e-07 ***
Neck        -0.014175   0.059800  -0.237   0.8129    
Abdomen      0.137760   0.022117   6.229 3.57e-09 ***
Hip         -0.002438   0.036103  -0.068   0.9462    
Thigh        0.145343   0.035583   4.085 6.79e-05 ***
Knee        -0.169290   0.071710  -2.361   0.0194 *  
Ankle        0.020266   0.066001   0.307   0.7592    
Biceps       0.042724   0.046566   0.917   0.3602    
Forearm      0.122550   0.061477   1.993   0.0478 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.043 on 170 degrees of freedom
Multiple R-squared:  0.8924,    Adjusted R-squared:  0.8867 
F-statistic: 156.6 on 9 and 170 DF,  p-value: < 2.2e-16

Możemy zauważyć, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę.

Code

summary(dane3$BMI-model3$fitted.values)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-11.81280  -3.07204   0.14458   0.06521   2.98281  10.48664 

Code

autoplot(model3)

Wykresy diagnostyczne wyglądają lepiej, niż w poprzednim modelu, spodziewamy się, że większość zalożeń będzie spełnionych.

Założenia Gaussa Markowa

Code

gauss_markov(model3)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.682	jednorodność wariancji
Test Goldfelda-Quandta	0.933	jednorodność wariancji
Test Durbina-Watsona	0.787	brak autokorelacji reszt
Test Breuscha-Godfreya	0.197	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.867	rozklad reszt jest normalny
Test Kolmogorova-Smirnova	0.98	rozklad reszt jest normalny
Test Lillieforsa	0.843	rozklad reszt jest normalny
Test RESET Ramseya	0.157	zależność jest liniowa
Test Rainbow Uttsa	0.858	zależność jest liniowa
Test Harveya-Colliera	0.908	zależność jest liniowa

W kolejnym kroku po zbudowaniu pełnego modelu przejrzymy się czy da się nasz model ulepszyć, widzimy, że model posiada 9 zmiennych, co jest sporo, chcemy zredukować je ilość za pomocą metody PCA, którą przedstawiamy poniżej:

Najpierw sprawdzamy korelację zmiennych:

Widać, że niektóre zmienne są mocno skorelowane, więc PCA będzie tutaj dobrym wyborem:

Code

model3.pca <- prcomp(uczacy[,-c(1:6,16:20)], scale = T)
summary(model3.pca)

Importance of components:
                          PC1     PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     2.5164 0.86709 0.7794 0.6490 0.53794 0.49801 0.43170
Proportion of Variance 0.7036 0.08354 0.0675 0.0468 0.03215 0.02756 0.02071
Cumulative Proportion  0.7036 0.78715 0.8547 0.9014 0.93360 0.96116 0.98187
                           PC8     PC9
Standard deviation     0.30138 0.26899
Proportion of Variance 0.01009 0.00804
Cumulative Proportion  0.99196 1.00000

Code

fviz_eig(model3.pca, addlabels = TRUE)

Code

fviz_contrib(model3.pca, choice = "var", axes = 9)

Widzimy, że pierwsze 3 składowych już dałyby nam niemal 86% wyjaśnionej wariancji, co byłoby zgodnie z kryterium.

Tworzymy zmienne PC1,PC2,PC3.

Code

PC1 <- model3.pca$x[,1]
PC2 <- model3.pca$x[,2]
PC3 <- model3.pca$x[,3]

Code

model3_2 <- lm(BMI~PC1+PC2+PC3, data=uczacy)
summary(model3_2)

Call:
lm(formula = BMI ~ PC1 + PC2 + PC3, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8746 -0.7642  0.0386  0.8813  4.0660 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 24.13333    0.08748 275.868  < 2e-16 ***
PC1         -1.09812    0.03486 -31.500  < 2e-16 ***
PC2         -0.80886    0.10117  -7.995 1.65e-13 ***
PC3         -0.42268    0.11255  -3.756 0.000235 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.174 on 176 degrees of freedom
Multiple R-squared:  0.8588,    Adjusted R-squared:  0.8564 
F-statistic: 356.8 on 3 and 176 DF,  p-value: < 2.2e-16

Code

gauss_markov(model3_2)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.355	jednorodność wariancji
Test Goldfelda-Quandta	0.949	jednorodność wariancji
Test Durbina-Watsona	0.872	brak autokorelacji reszt
Test Breuscha-Godfreya	0.371	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.334	rozklad reszt jest normalny
Test Kolmogorova-Smirnova	0.362	rozklad reszt jest normalny
Test Lillieforsa	0.433	rozklad reszt jest normalny
Test RESET Ramseya	0.03	brak zależności liniowej
Test Rainbow Uttsa	0.396	zależność jest liniowa
Test Harveya-Colliera	0.587	zależność jest liniowa

R^2 tego modelu to 0.85 Założenia Gaussa-Markowa z różnych testów troszkę się polepszyły.

Model zagnieżdżony

Code

mod_0 <- lm(BMI~1, data = uczacy)
met_tyl <- stats::step(model3, scope = c(mod_0, model3),
                  direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients

(Intercept)       Chest     Abdomen       Thigh        Knee     Forearm 
 -9.4910775   0.1444103   0.1339823   0.1549096  -0.1639500   0.1378829 

Code

model_zag <- lm(BMI~Chest+Abdomen+Thigh+Knee+Forearm, data=uczacy)
summary(model_zag)

Call:
lm(formula = BMI ~ Chest + Abdomen + Thigh + Knee + Forearm, 
    data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6718 -0.6431  0.0511  0.6350  3.5310 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.49108    1.51928  -6.247 3.11e-09 ***
Chest        0.14441    0.02415   5.981 1.23e-08 ***
Abdomen      0.13398    0.01978   6.774 1.85e-10 ***
Thigh        0.15491    0.02801   5.531 1.15e-07 ***
Knee        -0.16395    0.06147  -2.667  0.00837 ** 
Forearm      0.13788    0.05417   2.545  0.01178 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.033 on 174 degrees of freedom
Multiple R-squared:  0.8918,    Adjusted R-squared:  0.8887 
F-statistic: 286.7 on 5 and 174 DF,  p-value: < 2.2e-16

Z podsumowania można wnioskować, że prawdziwe BMI odchyla się od linii regresji średnio o około 1 jednostkę, tak samo jak w modelu pełnym.

Code

summary(uczacy$BMI-model_zag$fitted.values)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-2.67183 -0.64313  0.05105  0.00000  0.63502  3.53102 

Założenia Gaussa Markowa

Code

gauss_markov(model_zag)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.638	jednorodność wariancji
Test Goldfelda-Quandta	0.918	jednorodność wariancji
Test Durbina-Watsona	0.787	brak autokorelacji reszt
Test Breuscha-Godfreya	0.63	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.741	rozklad reszt jest normalny
Test Kolmogorova-Smirnova	0.979	rozklad reszt jest normalny
Test Lillieforsa	0.808	rozklad reszt jest normalny
Test RESET Ramseya	0.277	zależność jest liniowa
Test Rainbow Uttsa	0.813	zależność jest liniowa
Test Harveya-Colliera	0.945	zależność jest liniowa

Code

autoplot(model_zag)

Code

anova(model3, model_zag)

Analysis of Variance Table

Model 1: BMI ~ Chest + Neck + Abdomen + Hip + Thigh + Knee + Ankle + Biceps + 
    Forearm
Model 2: BMI ~ Chest + Abdomen + Thigh + Knee + Forearm
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    170 184.79                           
2    174 185.82 -4   -1.0382 0.2388 0.9161

Na podstawie analizy wariancji przyjmujemy, że model prostszy jest lepszy.

Predykcja

Code

pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,19)]

Code

cbind(pred, fakt = fakt$BMI) %>% 
ggplot(aes(x = pred, y = fakt))+
  geom_point()+
  geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Podsumowanie

Postać modelu:

\hat {BMI} = -9.49 + 0.14*Chest + 0.13*Abdomen + 0.15*Thigh - 0.16*Knee + 0.14*Forearm

Nasz model wyjaśnia niemal 89% wariancji. Zmienne wchodzące w skład modelu to obwody: klatka piersiowa, brzuch, uda, kolana i przedramię. Wskazuje to na dużą zależność tych miar ze wskaźnikiem masy ciała. Oznacza to, że dla przeciętnych mężczyzn wskaźnik ten jest wystarczająco dobrą miarą aby określić swój stan fizyczny. Pamiętajmy jednak, że przy dużej objętości tkanki mięśniowej nie jest on miarodajny i raczej wybralibyśmy procentowy wskaźnik poziomu tłuszczu.

Hipoteza 2

Czy poszczególne miary: obwody klatki piersiowej, brzucha, bioder, uda, kolana, kostki, bicepsa, przedramienia, nadgarstka, wiek, waga i wzrost mają wpływ na poziom tłuszczu?

Code

ggplot(dane3,aes(x=BodyFat,y=Weight2)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i wagi") +
  geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 250)

`geom_smooth()` using formula = 'y ~ x'

Code

ggplot(dane3,aes(x=BodyFat,y=Abdomen)) + geom_point() + theme_bw() + ggtitle("Zależność % tkanki tłuszczowej i obwodu brzucha") +
  geom_smooth(method='lm') + ggpubr::stat_regline_equation(label.x = 5, label.y = 115)

`geom_smooth()` using formula = 'y ~ x'

Model pierwszy

\hat{BodyFat} = \beta_{0} + \beta_{1} \cdot Chest + \beta_{2} \cdot Abdomen + \beta_{3} \cdot Hip + \beta_{4} \cdot Thigh + \beta_{5} \cdot Knee + \beta_{6} \cdot Ankle + \beta_{7} \cdot Biceps + \beta_{8} \cdot Forearm + \beta_{9} \cdot Age + \beta{10} \cdot Weight2 + \beta{11} \cdot Height2+\beta{12} \cdot Wrist

Code

model <- lm(BodyFat~Chest+Abdomen+Hip+Thigh+Knee+Ankle+Biceps+Forearm+Age+Weight2+Height2+Wrist, data=uczacy)
summary(model)

Call:
lm(formula = BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee + 
    Ankle + Biceps + Forearm + Age + Weight2 + Height2 + Wrist, 
    data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.1699 -3.0621 -0.0413  2.8888  9.0910 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -24.67529   25.91811  -0.952  0.34245    
Chest        -0.05061    0.11979  -0.422  0.67321    
Abdomen       0.81385    0.10564   7.704 1.11e-12 ***
Hip          -0.10127    0.15618  -0.648  0.51761    
Thigh         0.32767    0.16133   2.031  0.04384 *  
Knee          0.32390    0.31165   1.039  0.30017    
Ankle        -0.10546    0.28083  -0.376  0.70774    
Biceps        0.21557    0.18828   1.145  0.25388    
Forearm       0.37294    0.25023   1.490  0.13801    
Age           0.06188    0.03637   1.701  0.09073 .  
Weight2      -0.21411    0.17333  -1.235  0.21847    
Height2      -0.07795    0.08480  -0.919  0.35931    
Wrist        -1.95932    0.61351  -3.194  0.00168 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.097 on 167 degrees of freedom
Multiple R-squared:  0.7645,    Adjusted R-squared:  0.7476 
F-statistic: 45.17 on 12 and 167 DF,  p-value: < 2.2e-16

Code

summary(uczacy$BodyFat-model$fitted.values)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-9.16994 -3.06208 -0.04125  0.00000  2.88881  9.09097 

Założenia Gaussa Markowa

Code

gauss_markov(model)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.866	jednorodność wariancji
Test Goldfelda-Quandta	0.578	jednorodność wariancji
Test Durbina-Watsona	0.987	brak autokorelacji reszt
Test Breuscha-Godfreya	0.859	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.226	rozklad reszt jest normalny
Test Kolmogorova-Smirnova	0	rozklad reszt nie jest normalny
Test Lillieforsa	0.359	rozklad reszt jest normalny
Test RESET Ramseya	0.352	zależność jest liniowa
Test Rainbow Uttsa	0.845	zależność jest liniowa
Test Harveya-Colliera	0.641	zależność jest liniowa

Model zagniezdżony

Code

mod_0 <- lm(BodyFat~1, data = uczacy)
met_tyl <- stats::step(model, scope = c(mod_0, model),
                  direction = 'backward', test = 'F', trace = 0)
met_tyl$coefficients

 (Intercept)      Abdomen        Thigh      Forearm          Age      Weight2 
-44.07772256   0.83167546   0.41690583   0.51621459   0.07463536  -0.30133107 
       Wrist 
 -1.82469765 

Code

model_zag <- lm(BodyFat~Abdomen+Thigh+Weight2+Wrist, data=uczacy)
summary(model_zag)

Call:
lm(formula = BodyFat ~ Abdomen + Thigh + Weight2 + Wrist, data = uczacy)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.5177 -2.8102 -0.0278  3.1434  9.9275 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -41.85380    9.78803  -4.276 3.12e-05 ***
Abdomen       0.91504    0.06283  14.563  < 2e-16 ***
Thigh         0.34524    0.12056   2.864 0.004701 ** 
Weight2      -0.31348    0.08241  -3.804 0.000196 ***
Wrist        -1.08634    0.52264  -2.079 0.039118 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.143 on 175 degrees of freedom
Multiple R-squared:  0.7476,    Adjusted R-squared:  0.7419 
F-statistic: 129.6 on 4 and 175 DF,  p-value: < 2.2e-16

Możemy zauważyć, że prawdziwy poziom tłuszczu odchyla się od linii regresji średnio o około 4 punkty procentowe.

Code

summary(uczacy$BodyFat-model_zag$fitted.values)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-9.51767 -2.81022 -0.02781  0.00000  3.14340  9.92750 

Code

autoplot(model_zag)

Założenia Gaussa Markowa

Code

gauss_markov(model_zag)

Założenia Gaussa Markowa

Nazwa testu	P-wartość	Wniosek
Test Breuscha-Pagana	0.811	jednorodność wariancji
Test Goldfelda-Quandta	0.455	jednorodność wariancji
Test Durbina-Watsona	0.974	brak autokorelacji reszt
Test Breuscha-Godfreya	0.168	brak autokorelacji do rzędu 3
Test Shapiro-Wilka	0.334	rozklad reszt jest normalny
Test Kolmogorova-Smirnova	0	rozklad reszt nie jest normalny
Test Lillieforsa	0.25	rozklad reszt jest normalny
Test RESET Ramseya	0.052	zależność jest liniowa
Test Rainbow Uttsa	0.885	zależność jest liniowa
Test Harveya-Colliera	0.718	zależność jest liniowa

Zauważamy, że R dopasowane jest zbliżone w modelu zagnieżdżonym do R dopasowanego w modelu pełnym. Oba modele spełniają wszystkie założenia, jedynie test Kołmogorowa - Smirnova wykazuje brak normalności rozkładu reszt. Wszystkie zmienne wchodzące w skład modelu zagnieżdżonego są istotne statystycznie, więc skłanialibyśmy się ku modelowi zagnieżdżonemu.

Code

anova(model, model_zag)

Analysis of Variance Table

Model 1: BodyFat ~ Chest + Abdomen + Hip + Thigh + Knee + Ankle + Biceps + 
    Forearm + Age + Weight2 + Height2 + Wrist
Model 2: BodyFat ~ Abdomen + Thigh + Weight2 + Wrist
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    167 2803.8                           
2    175 3004.4 -8   -200.65 1.4939 0.1628

Na podstawie analizy wariancji stwierdzamy, że model prostszy jest lepszy.

Predykcja

Code

pred <- round(predict(model_zag,testowy),1)
fakt <- testowy[,c(1,3)]

Code

cbind(pred, fakt = fakt$BodyFat) %>% 
ggplot(aes(x = pred, y = fakt))+
  geom_point()+
  geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Podsumowanie

Postać modelu:

\hat {BodyFat} = -41.85 + 0.92*Abdomen + 0.35*Thigh - 0.31*Weight2 - 1.09*Wrist

Nasz model wyjśnia ponad 74% wariancji. Miary wchodzące w skład modelu to obwody: brzuch, uda, nadgarstek oraz waga w kilogramach, dzięki czemu każdy (mężczyzna) może łatwo sobie wyznaczyć procentowy poziom tłuszczu i sprawdzić do jakiej kategori należy. Pomaga to w kontrolowaniu kondycji fizycznej, która ma wysoki wpływ na zdrowie człowieka.